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Abstract 

The factorizable discretization of Sidilkovei for 
the compressible Euler equations previously demon- 
strated for channel flows has been extended to ex- 
ternal flows . The dissipation of the original scheme 
has been modified to maintain stability for mod- 
erately stretched giids. The discrete equations are 
solved by symmetric collective Gauss-Seidel relax- 
ation and FAS multigrid. I nlike the earlier work, 
ordering ; the grid vertices in the flow direc tion has 
been found to be unnec cssary. Solutions for essen- 
tial incompressible flow (Mach 0.01) and supercrit- 
ical flows have been obtained for a I\ arm an- 1 refltz 
airfoil with a conformally mapped grid, as well as 
a A A ( W 0012 on an algebraically generated grid. 
The current work demonstrates nearly C)( n ) con- 
vergence for subsonic and slightly transonic flows. 

Introduction 

In the past several years, there has been an ongo- 
ing effort in the Computational Modeling and Simula- 
tion Branch of NASA Langley Research Center, in col- 
laboration with researchers at H ASH, to develop algo- 
rithms demonstrating "Textbook Multigrid Efficiency" 
(TME)*’ The central theme of this effort lias been 
the idea of Brandt ^ that ideal convergence rates for the 
numerical solution of part ial differential equations can be 
achieved if the relaxation scheme distinguishes between 
elliptic, parabolic, and hyperbolic partitions of the dif- 
ferential operator. Construction of an ideally converging 
method requires that each of these partitions be treated 
independently and optimally. To achieve this, efforts are 
being made to construct so-called "fact orizable" differ- 
ence schemes for the equations of fluid dynamics. A fac- 
torizable scheme is one in which the discretization cor- 
rectly separates the elliptic, parabolic, and hyperbolic 
partitions. 

In t he case of steady inviscid flow, the Euler equations 
can be considered as a composition of two subsystems. 
One subsystem corresponds to the equations governing 
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entropy and vorticity advert ion. This subsystem is hy- 
perbolic in space. The ot her subsystem corresponds to a 
full potential operator, which is elliptic for subsonic flow 
and hyperbolic lor supersonic flow. For a purely super- 
sonic flow, space marching is the most efficient way of 
solving the Euler equations. For subsonic flow, the el- 
lipticity of the full potential factor should be effectively 
handled by multigrid. However multigrid is not effective 
for advert ion operators, as the coarse grid only gives part 
of the correction for certain smooth components of the 
error. Existing multigrid methods for subsonic and tran- 
sonic flow rely on the coarse grid to smooth the entire 
system. As such, they are fundamentally limited by the 
ineffectiveness of the coarse grid in correcting the part ol 
the error corresponding to advection factor. This same 
difficulty is true for high Reynolds number viscous flows. 

Sidilkovei J has devised a ('artesian grid discretiza- 
tion of the Euler equations that separates these* two 
subsystems. In a previous work, Roberts, Sidilkovei 
and Thomas'* extended this scheme in conservation form 
to general body-fitted curvilinear coordinates. Solutions 
to the equations were obtained using collective point 
Gauss-Seidel relaxation as a smoother combined with 
FAS multigrid. Results for subsonic and transonic chan- 
nel flows were presented in Ref. 3. which demonstrated 
essentially grid independent convergence rates. All the 
grids used for these cases were quasi-uniform, with grid 
cells of unit aspect ratio. Furthermore, the channel wall 
had a continuous curvature to preclude the presence of 
stagnation points, anti the relaxation sweep was made in 
the flow direction. 

in t he present paper, the work of Roberts, et al. is ex- 
tended to external lifting flows around airfoils with bot h 
subcritical and supercritical freestream Mach numbers. 
The presence* of a stagnation point at the airfoil lead- 
ing edge requires an additional dissipation term in lilt* 
momentum equation. The form of this dissipation was 
presented in Ref. 3 and is found to be effective. Compu- 
tations are shown in which an odd-even vort icity error is 
propagated downstream of the airfoil. This error origi- 
nates at the trailing edge of the airfoil, and causes erratic 
convergence. A dissipation term to smooth this vortical 
error is added to the momentum equations. 

The mathematical formulation of the scheme, which 
has been presented by Sidilkovei 0 and Roberts, et al.’* is 
briefly recounted in the following section. Modifications 
to the scheme as presented in Ref. 3 are noted. Follow- 



ing that, the multigrid relaxation is described. Subsonic 
channel flow results are presented to demonstrate the 
second-order convergence of the discretization. Results 
for subsonic and transonic airfoils are present ed next . 
Subsonic solutions for a Karinan-Trefftz airfoil are used 
to compare with an exact solution. In addition, a tran- 
sonic solution for this airfoil is presented. The Karinan- 
Trefftz airfoil airfoil uses a conformally-mapped grid with 
unit aspect ratio ceils. To demonstrate the effectiveness 
of the scheme when applied to a more realistic case, sub- 
sonic ami transonic flow around a NA(-A 0012 is also 
presented. These computations are performed on a grid 
generated algebraically, with moderate grid stretching, 
skewness, and aspect ratio. 

Mathematical Formulation 

The artificial dissipation of the factorizable scheme 
can be described by first presenting the modified equa- 
tion, or first differential approximation (FDA), of the 
discrete scheme. This is the differential equation which 
is found by expanding the difference equation in a Taylor 
series about each grid vertex and considering the lead- 
ing terms of the truncation error. These terms are t he 
artificial dissipation of the scheme. 

The starting point, for the scheme is the two- 
dimensional Euler equations in non-conservation form. 
Let p be the density, u = w 4- jv be the velocity, and p 
be the pressure. The entropy is defined as 



where po and po are a reference pressure and density, 
respectively, and 7 is the ratio of the specific heats. 
T hen the Euler equations may be written in the vari- 
ables (.s, u,c,p): 


it -Vs = 0 , 

(2a) 

u-Vti + - vp = o. 

(2b) 

p 


pc 2 V- i 7 + u V p — 0. 

(2c) 


The factorizahiiitv of the scheme depends on the form 
of the artificial dissipation added to this system of equa- 
tions. 

The entropy is only weakly coupled to the momen- 
tum and the pressure equations through the equation 
of state (1). In fact, the entropy equation corresponds 
to one of the advection factors of the Euler equations. 
Therefore, it may be discretized independently of the 
momentum and pressure equations in any appropriate 
way without affecting the factorizability of the scheme. 
The advection operator u V uses simple upwind differ- 
encing in Eq. (2a). Let (£,?/) be a general curvilinear 
coordinate system and define the contravariant compo- 
nents of the velocity, (b\ V) by the transformation 



In this coordinate system TV ?7 = l c)^ + Vd, t . The equa- 
tions are discretized on a uniform grid in (£,?/) space, 
with a grid spacing A£ = A/7 = 1. The FDA of the 
first-order upwind difference approximation to w-V is 

<1 = w-V-i | / | ^ H’l^ = (), (4) 

and the entropy equation is discretized as 

q* = 0. ( r >) 

A second-order upwind discretization of the advection 
operator has also been used in Eq. ( 5 ). However, it was 
shown in Ref. 3 that the use of a second-order advection 
operator has an insignificant affect on the computed re- 
sults. 

The dissipation for the momentum and pressure equa- 
tions, Eqs. (2b) and (2c) given in Ref. 3 is the multidi- 
mensional upwind dissipation of Sidilkover^. In vector 
notation, this dissipation is written as 

.7Vt7 + -Yp - —V (pc 2 Yu + it-Yp) = 0. (6) 

P -P ( 

pr‘Y- it -j- (7 V }> — pc— ^(7-Vi7-f — = 0, (7) 

where c is the speed of sound, <r„, and are scaling 
coefficients, and ( is a lengt h scale proportional to the 
grid spacing. Note that the dissipation of the momen- 
tum equation is the gradient of the pressure equation 
residual, and the dissipation of the pressure equation is 
the divergence of the momentum equat ion residual. Also 
note that the curl of Eq. (6) is identical to the curl of 
Eq. (2b), i.e., vorticity equation of the governing Euler 
equations is unaffected by the artificial dissipation. The 
vorticity equation corresponds to the second advection 
factor of the Euler equations. 

If the advection operator, pressure gradient, and di- 
vergence term in Eqs. (G) and ( 7 ) are discretized us- 
ing central differences, the scheme is second-order accu- 
rate and factorizable. However, such a scheme is not 
/^-elliptic. This lack of //-eilipticity is a result of the 
central difference approximation to the advection term 
in the momentum equation, Eq. (6). Sidilkover shows 
that this advection operator corresponds to vorticity 
advection 0 . Replacing this operator by the first-order 
upwind approximation in Eq. ( 4 ) restores /j-ellipticity 
and the scheme remains factorizable, but it now becomes 
only first-order accurate. Note that the pressure advec- 
tion term in Eq. ( 7 ) continues to be approximated by 
second-order central differences. 

To obtain second-order accuracy while maintaining 
h-ellipticity, the advection operator in the momentum 
equation must be replaced with an upwind difference op- 
erator. However, this operator also contributes to th< full 
potential operator. Simply replacing the central differ- 
ence with an upwind difference will change the form of 
the full potential operator, which we wish to avoid. To 
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obtain a second-order upwind advert ion operator with- 
out changing the full potential operator, Eq. (0) is mod- 
ified to be 


u Vit + — Vp H u xVx(a V u H — V 7 p 

p - ( i V p 

- (pr a V-ff + u-Yp) = 0. (8) 

where q = niax( |u|, |r| ). Because the additional term 
involves the curl of tlie momentum equation, it is trans- 
parent to irrotational disturbances. If a solution is a 
pure potential solution, then not only is Eq. (G) satisfied 
identically, but Eq. (K) is satisfied as well. Because the 
continuity equation (7) is unchanged, a potential solu- 
tion to Eq. ((>) and (7) is a solution to Eq. (X) and (7). 
In other words, the modificat ion to the momentum equa- 
tion does not change the full potential operator, which 
is what is wanted. 

Expanding the additional term in Eq. (X) the principal 
part of the operator is 

— u x V x (it -V it H — Vp 

iq V p 

j- L(«.V) 2 h + V ('<7 V “ + - ^ +Sxvbvp . 

(*>) 

The first term in the square brackets may be combined 
with the advert-ion operator in Eq. < K ) to yield the 
FDA of a second-order upwind advert ion operat or. This 
second-order operator operator is denoted as quo- The 
pressure terms can be manipulated such that the remain- 
ing terms in the square brackets be can be written in the 
form V/). This puts the momentum equation in the form 

qn„M + Yl) + -Yp - — Y (pr 2 V ii+u Vp) = 0. 

p >pr 

( 10 ) 

In Roberts, et al4 the pressure term in E<|. (9) was 
ignored because the principal part of this term vanishes. 
In fact, if the flow is barotropic the term vanishes iden- 
tically. However, failure to include this term in the con- 
servative discretization gives rise to a first order error, a 
point that was not recognized at the time. In the current 
work these terms are included in the discretization. 

Discretization 

The factorizability of the FDA is a necessary condi- 
tion for the factorizability of the difference 1 scheme, but 
it is not sufficient. For the difference scheme to be fac- 
torizable. the difference operators must commute in the 
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same way as the differential operators' . Introducing the 
difference approximations to the partial derivative oper- 
ators. 



Dc — ()$ + • • • , c) f} — <) t f 4- ■ • ■ , 

() f ^c = c)i 4 * • • • , — <K t + • • • « 
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the following conditions must hold 

Jt Jt Jt Jt 

f ss ( V/ “ 

■ Jt . Ji ■ Jt Jt 

= (> (’A ■ 

■ Jt Jt .■Jt 


The following difference operators satisfy this condition: 
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To write the complete discrete scheme, the subscript h 
is used to denote 1 a standard difference to tl«* corre- 
sponding operator, anti the addition of the overbar ( ) 
denotes the “wide" differences given above*. The* second- 
difference expressions may be expressed in flux form by 
taking a six-point difference centered on an edge between 
two vertices, and then taking a two-point difference of 
those expressions centered on a vertex. The* subscripts t 
and v are used to denote difference operators centered 
on an edge or a vertex, respectively. The fully discrete 
scheme is then written as 


Jt 

l h’, - 


h M 

<1 s = 0. 


(Ha) 


qn t ,u + V,l) h + -V h p 
P 

- V* (pJv’MT + u-Y'.’p'j =0. (II h) 
pr 2 V h -u + uY'’p 

- pr^-Y'l- (it-Yl'u + J- (y 1 : + t 7 ?) p'j =0. (I Ir) 

Details of t he velocity terms in D are found in t he Ref. 3. 
The pressure terms in f) an* evaluated using the with* 
difference operators defined above*. 

The scaling coefficients are 

= max( M , M r ), a f , = y - V . (12) 

max( .1/, a/, ) 

where M is the local Mach number, and M, is a cutoff 
Mach number to prevent division by zero, The cutoff is 
chosen to be (){h), and essentially becomes active near 
stagnation ]>oints. The pur]>ose of the rescaling of the 
pressure equation dissipat ion. cr v > is lo prevent the ellip- 
tic factor in the discrete equations from becoming the 
skewed Laplacian operator in tire incompressible limit. 



Currently, we take M, — |AA/|, where v — 2 

and AM is the two-point difference in Mach nuniher 
on an edge. When M r > M. it is necessary to add ad- 
ditional dissipation to the advection operator q h . The 
form of this dissipation is a five-point pseudo- Laplacian, 

d.* v it = — max (0, M, — M) ~j + 0? } ) w, (13) 

wliere J is the Jacobian of the coordinate transforma- 
tion, J = j zfjjj — UzJ'ii- This operator is added to botli 
the entropy and the moment um equation, and is cast in 
flux form in order to maintain conservation. No attempt 
has been made to optimize either the coefficient u or the 
form of the operator <T P . 

When solving for the flow around an airfoil, it has 
been observed that a short -wavelength vorticity error is 
generated near the trailing edge and is propagated down- 
stream. This error is an odd-even error of the stream wise 
velocity component in the crossflow direction. The error 
persists to the outflow boundary, and can cause conver- 
gence difficulties. An additional dissipative term must 
be added to the momentum equation to reduce the er- 
ror. ( are must be taken to insure that the dissipation 
does not affect the fact orizabi lily of the scheme and pre- 
serves second- order accuracy. 

The form of this dissipation follows from a formulation 
of Giles' . who observed that V x £ — V(V-t7) — V 2 u. 
The Laplacian acting on u is dissipative. In the present 
scheme, the addition of a term of this form to the mo- 
mentum equation will not affect the factor izabilitv of the 
scheme, as it is transparent to the irrotational, potential 
part of the solution. This vorticity dissipation is 

d,.n = /< | a- x Au| V' 1 x . ( 14) 

where p is a dissipation coefficient, An is ati undivided 
difference of the velocity in a grid direction, ( is the grid 
spacing directed in that direc tion, and Jv is the discrete 
vorticity centered on an edge* and computed using the 
wide differences. Currently the coefficient p is taken to 
be 0.1. No attempt has been made to optimize this dis- 
sipation. 

As written, the Eq. (11) is valid only for subsonic 
flows. This is because the pressure differences in the 
artificial dissipation terms are not fully upwinded in a 
supersonic zone. A simple modification to Eq. (11) can 
be made by rescaling the gradients of the pressure when 
the flow becomes sonic. Introducing the parameter k, 
defined as 

k = max( 1 , A/ 2 ), (1 5 ) 

t he final form of the scheme is 

q /} .s — a = 0, (16a) 

q hu « - d*pu + d v u + V l) h -f ^ Y h p 

_ l 2 v;-.,7 + = o, (i6b» 

2 pc V K 


2 T— 7 h -» , -* t — r h 

pc V • tt + u • V p 

- pc^v 1 !- (u Vi’ii + (v 1 : + VH) ,)j = 0. ( 16c) 

Because the difference equations (16a), (16b) 

and (16c). can be written as a central difference part 
plus a dissipation, it is straightforward to obtain a con- 
servative discretization. The conservation form of the 
Euler equations are discretized using a cent ral-difference 
finite- volume approximation. 


Y h ‘(pu) = (), 

(17a 

,h (puu) + Y h p = 0. 

(]"l> 

+ p)«i] = o 

(!"<• 


where t = u-u/2 -f p/{p(l - 1)) is the total energy. 
The dissipative fluxes on each cell face are computed 
in terms of (s, u,e,p) using the appropriate difference 
operators in Eq. (16). The artificial dissipation in 
Eqs. (16b) and (16c) can be rewritten as 

A <7 = — V r L) h + dgpii — d v tt 

+ vj! [— frV'/ u + — «V?A1 . (IS) 

S t> = V*. [pru-V h n + ™ (y 1 : + VJ ) p)] . ( 19) 

This is a conservative form of the dissipation cast in 
terms of the primitive variables. The terms inside the 
square bracket in Eqs. ( 18) and (19) are now interpreted 
as dissipative fluxes. These terms are discretized on cell 
faces according to the appropriate wide or narrow dif- 
ferences. The first and second-order upwind advection 
operators t\ h and qj^, and the terms V v D h in Eq. (18) 
may be recast in conservation form and evaluated on 
cell faces. The length scale f is evaluated on the cell face 
and is taken as either or k 7] is used, whichever is more 
nearly aligned with the flow direction. This assures that 
the advection operator in Eq. (16b) is fully upwinded 
when the flow is grid-aligned. The scaling coefficients a m 
and <t p are evaluated using the Mach number on the face. 

In addition to the dissipation terms, examination of 
Eqs. (1Gb) and (16c) shows that the pressure gradi- 
ent terms are discretized using the wide stencils. Tiie 
central-difference part of the conservation equations (17) 
must be corrected to account for these wide differences. 
This is done by adding a term to the dissipative fluxes 
for the momentum and pressure equations as follows: 

SCu-Hu- - (V* - V h }p, (20) 

5p <$p — w- (v A — V /! ) p. (21) 

Once the dissipative fluxes (A.s, An, Ar. Ap) have been eval- 
uated, they are converted to the conservation variables 
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+ p/p is the total out 

tialpy. 




Solution Procedure? 

The solution of the discrete equations is performed us- 
ing a symmetric point, collective Gauss-Seidel iteration, 
which lias been found to be a very effective smoother. 
In the previous work‘d the grid vertices were ordered in 
the flow direction, so that the advection factor would be 
solved by marching during the forward sweep. Similarly, 
in the present work the grid vertices are sorted in the 
freest ream direction. Each Gauss-Seidel sweep is under- 
relaxed with a factor of 0.7 for stability. The residuals 
of the conservation equations ( r /M r,„ , r,„ ) are com- 
puted at each vertex using the discretization presented 
in the previous section. These residuals are then trans- 
formed to the residuals of the primitive variables by the 
inverse of the transformation in Eq. (22). 

At each point, the update to the solution is given by 


/a.A 


f r \ 


= - 

r n 

r, 

W 

> J 

vJ 


where M is the matrix of coefficients of the primitive 
variables at vertex The coefficients are found by 

collect ing the contributions to vertex (/, j) from the dissi- 
pation terms on t he four surrounding faces. Because the 
entropy equation decouples from the rest, of the system, 
t his is a block diagonal matrix where the upper block is 
the upper left-hand entry, and the lower block is a 3 x 3 
mat rix of coefficients multiplying ih.j, v tJ , an dp,j. This 
matrix is easily inverted. 

The relaxation is accelerated using a standard Full- 
Approximation Scheme (FAS) multigrid. A sequence of 

grids Gj v , (7/v-i Go is used, where G/ v is the finest 

and Go the coarsest. Let L a _i be the coarse grid oper- 
ator. u be the vector of the conservation variables, /*_ | 
be the fme-to-coarse grid rest riction operator, and /£ _1 
be the coarse- to- fine grid prolongation operator. If u* is 
the current solution on grid k\ the residual on this grid 
is i*a = f A — L a u a-. This is the residual of the conserva- 
tion equations, not the primitive equations. This leads 
to the coarse- grid equation 

La_iIU-i = f*A — l = / A _i I'a -b La-i j • (23) 

After solving the coarse-grid equation for ua_i, the fine- 
grit l solution is corrected by 

t*A *U -f- 1 — /Ldu) - (25) 



1 igure 1. ( 'divergence of lift, drag, and entropy for 
channel flow with a 10 ( /{ bump on the lower wall 


Equal ion (23) is solved by applying the same relaxation 
procedure that is used to solve the fine-grid equation. 
Multigrid is applied recursively to the coarse-grid equa- 
tion. On the coarsest grid, many relaxation sweeps are 
performed to insure that the equation is solved com- 
pletely. A conventional \ cycle or U -cycle is used. 

Results 

To confirm the order of accuracy of the scheme, solu- 
tions for the channel configuration used in Ref. 3 were 
obtained for a sequence of grids starting at 7 x 3 and 
ending at 769 x 257. Results are shown in Fig. 1. The 
lilt and drag values are the integrated forces on the lower 
wall. Richardson extrapolation was used to compute the 
exact lift coefficient. The Li norm of the difference be- 
tween t he local and freest ream ent ropy is also shown. 
It is seen that all three quantities exhibit second-order 
convergence in the grid spacing. 

The first set of external flow results are for flow around 
a symmetric Karman-Trefftz airfoil with a thickness ratio 
of about 0.15. An O-grid was generated using a confor- 
mal mapping and is shown in Fig. 2. The radial stretch- 
ing factor was chosen to maintain unit ceil aspect ratio, 
and grid lines are orthogonal. The finest grid for all cast s 
is 257 x 257 and the coarsest grid is 5 x 5. 'The outer 
boundary is approximately 150 chord lengths from the 
airfoil. 

Dirichict conditions are used at the outer boundary, 
a uniform freest ream being specified. No far- field vortex 
correction is applied. One-sided first-order differences 
are used at the wall. The tangency condition u • n = 0 is 
enforced explicitly by zeroing-out the component of the 
momentum equation residual in the direction normal to 
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Figure 2. Conformal grid for Karman-TrefFtz airfoil. 


the wall. 

Solutions are obtained using a FMG cycle to initial- 
ize the solution. A solution is computed on the coarsest 
grid and prolonged to the next grid to obtain a starting 
solution for that grid. This procedure is continued re- 
cursively until the finest, grid is reached. On each grid, 
a \ ( 1, 1 ) inultigrid cycle is used to solve the system. Ten 
mult igrid cycles are run on each of the coarse grids, and 
t wenty cycles are run on the finest grid. After each sym- 
metric Gauss- Seidel relaxation sweep, an additional six 
st ream wise sweeps are done on the wall and its neigh- 
boring row of vertices. 

f or the first case the freestream Mach number is 0.01 
and the angle-of-at tack is 2°. A comparison of the com- 
puted surface pressure coefficient with the exact incom- 
pressible solution is shown in Fig. 3. Excellent agree- 
ment is seen. Convergence rates for the continuity equa- 
tion residual as well as the lift and drag coefficients are 
shown in Fig. 4. On the coarse grids up to the 129 X 129 
grid the convergence rate is essentially independent ol 
the grid spacing, with an asymptotic rate of about 0.25. 
This slows down considerably on the 257 x'257 grid which 
has an asymptotic rate for the continuity residual of ap- 
proximately 0.G per cycle. There is some slightly erratic 
behavior in the convergence rate on the finest grid. Both 
lift and drag are converged after only two to four multi- 
grid cycles on each grid. 

A second case for the Karman-TrefFtz airfoil at the 
same freestream Mach number but an angle- of -at tack 
of 40° is shown in Figs 5 and 6. The comparison of sur- 
face pressures with exact values shows very good agree- 
ment. The 65 x 65 grid solution misses the suction peak 
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Figure 3. Surface* C p comparisons for a Karman- 
1 refit z airfoil, = 0.01, n = 2°, finest grid 257 x 

257 grid. 



Multigrid cycles 


Figure 4. Convergence of residual and forces for a 
Karman-TrelTtz airfoil, A/ irv = 0.01, o = 2°, finest 
grid 257 x 257 grid. 
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Figure 5. Surface C v comparisons for a Karman- 
Trefftz airfoil, — 0.01, n = 10°. finest grid 257 x 
257 grid. 

noticeably, but the two finest grid are very close to the 
analytic solution. The convergence rates shown in Fig. 6 
are very similar to the previous case. The asymptotic 
rate for the continuity equation residual is about 0.35 
per cycle on the coarser grids, and about 0.6 on the 
finest grid. The lift is converging more slowly on the 
finest grid, taking about five multigrid cycle to reach its 
asymptotic value. 

A transonic case is shown next, with a freest ream 
Mach number of 0.7 and an angle-of-attack of 1°. The 
surface C v distribution is shown in Fig. 7. showing the 
upper-surface shock at approximately 30% chord. The 
peak Mach number on the airfoil is approximately 1.15. 
Convergence rates are shown in Fig. X and are compa- 
rable to the subsonic case. The forces are converging in 
about three to four cycles on each grid. 

The grid for the cases shown above is ort hogonal with 
unit aspect ratio cells. To ascertain how the scheme be- 
haves in a more realistic case, subsonic and transonic 
solutions for a NACA 0012 have been obtained on an 
algebraically generated O-grid. 'The grid generator uses 
the transfmite interpolation method of Eriksson . The 
grid near the airfoil is shown in Fug. 9. The fine-grid 
dimensions are 257 x 129 with the outer boundary lo- 
cated approximately 100 chords from the airfoil surface . 
The grid is more skewed at the trailing edge than the 
Karman-Trefftz grid, and the cell aspect ratios vary from 
about 1:1 to 4:1. This grid is representative of a good- 
quality O-grid for a single-element airfoil. The coarsest 
grid is 9 x 5, and the same FM( ! cycle and boundary 
conditions are used as for the Karman-Trefftz airfoil. 

A solution for frees tream Mach number of 0.01 and 



Multigrid cycles 

Figure (>. (-divergence of residual and forces for a 
Karman-Trefftz airfoil. 3/^ = 0.01, n = 10°. finest 
grid 257 x 257 grid. 
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Figure 7. Surface ( f p comparisons for a Karman- 
Trefftz airfoil, = 0.70, n = 1°, finest grid 257 x 
257 grid. 
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Figure 8. Convergence of residual and forces for a 
Kami an-' Trefftz airfoil. M - v = 0.70. tv = 1°, finest 
grid 257 x 257 grid. 

angle-of-attack of 2° is shown in Fig. 10. < Convergence 
rates are shown in Fig. 11. A deterioration in the con- 
vergence rates compared to the Karman-Trefftz airfoil is 
observed, and the rate now appears to be slightly more 
grid-dependent . This is attributed to the varying aspect 
ratio of the grid cells. The asymptotic rate on the finest 
grid is about 0.5 per cycle. I lie forces are converging in 
about six to seven cycles on all grids. 

The final results is for the standard test case of AI iX . — 

0.80, a = 1.25°. This is a much more severe transonic 
case than the one shown for the Karman-Trefftz airfoil. 
The solution in Fig. 12 shows the strong suction-surface 
shock at about 60% chord, and the weak pressure- surface 
shock. The peak Mach number on the suction surface 
is about 1.35. The convergence rates for this case are 
shown in P ig. 13. The asymptotic rate for the continuity 
residual is quite slow on the finest grid, although the 
forces are converging in six to seven cycles. 

A major part of the slowdown in the convergence rate 
stems from the anisotropy of the full-potent ial factor in 
the transonic regime. This was not as significant as for 
t he Karman-Trefftz result shown in Fig. 8 because of the 
relatively small size of the supersonic region in that case 
and the lower freest ream Mach number. Point relaxation 
becomes decreasingly effective as a smoother for strongly 
supercritical flows. This can be alleviated by using line 
relaxation in the radial direction, which has the addi- 
tional benefit of being more effective on stretched grids. 

Conclusions 

A fact oriz able discretization of the compressible Eu- 
ler equations on general curvilinear body-fitted grids has 
been presented. This work is an extension of the previ- 
ous results of Roberts. Sidilkover and Thomas to lift - 



Figure 9. Algebraic grid for NA( A 0012 airfoil. 


ing airfoil flows. A modification of the original dissi- 
pation lias been introduced to preserve the stability of 
the point collective Gauss- Seidel smoother without de- 
stroying the factorizabiiity of the discretization A sig- 
nificant simplification of the original scheme has been 
the elimination of stream wise relaxation and replacing 
it with simple lexicographic symmetric relaxation. Both 
subsonic and transonic results for Karman-Trefftz and 
NACA 0012 airfoils demonstrates good accuracy and 
fast convergence. Nearly O(n) convergence rates for the 
residual are demonstrated subsonically and for slightly 
supercritical flows. The convergence for a more strongly 
transonic flow is significantly slower primarily because 
of the anisotropy of the differential operator in the tran- 
sonic regime. The use of line relaxation should improve 
this, and additionally is should be well-suited for more 
severe grid stretching. 
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